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FLAME  - A SLOW-FLOW  COMBUSTION  MODEL 


I.  Introduction 

This  paper  describes  FLAME,  a detailed,  two-dimensional,  time-dependent  numerical 
model  of  slow  flow  combustion  systems.  The  model  has  been  built  around  the  solution  of 
self-consistent  hydrodynamics  equations  in  the  ’slow  flow”  approximation  and  includes  the  basic 
transport  processes  of  molecular  diffusion  and  thermal  conductivity.  Other  physical  mechan- 
isms, such  as  detailed  chemical  reaction  schemes,  turbulence  models  and  heterogeneous 
phenomena  etc.;  can  be  included  on  a modular  basis.  This  facilitates  modifying  geometries  and 
modelling  various  processes  peculiar  to  specific  systems.  In  addition,  as  new  models  for  tur- 
bulence and  chemical  kinetics  become  available,  they  can  be  added  in  a relatively  straightfor- 
ward manner. 


The  basic  geometry  of  the  FLAME  model  is  cylindrical^  symmetric  (r  - z)'.  Other 
geometries  are  possible,  but  for  many  problems  of  interest,  namely  gas  jets,  diffusion  flames, 
fire  spread,  and  turbulent  mixing  studies,  this  would  appear  to  be  the  best  geometry.  The 
numerical  model  is  fully  two-dimensional  in  order  to  treat  quantitatively  processes  such  as 
buoyancy,  convection,  shear  flow  and  turbulence,  which  cannot  be  treated  sensibly  in  a one- 
dimensional calculation.  A fully  three-dimensional  calculation  using  the  same  techniques  is 
possible  but  would  be  significantly  more  expensive.  An  important  aspect  of  our  approach  has 
been  to  improve  the  resolution  of  sharp  gradients  via  Flux-Corrected  Transport.  Nevertheless, 
the  basic  model  is  Euierian  and  hence  resolution  of  gradients  is  necessarily  limited. 

The  set  of  equations  modelled  in  FLAME  is 
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This  is  a very  general  formulation.  The  approximations  and  assumptions  which  we  use  in  solv- 
ing these  equations  will  be  discussed  in  Section  III.  The  symbols  are 
p — total  mass  density  (gm/cm3) 
n - total  (particle)  density  (#/cm3) 
p,  — mass  density  of  the  ith  species  (gm/cm3) 
n,  — particle  density  of  the  ith  species  (#/cm3) 

V — bulk  fluid  velocity  (cm/sec) 

£ — V x V - vorticity  (sec-1) 

D — V ■ V - divergence  of  the  velocity  (sec-1) 

g — gravity  (-980  cm /sec2)  - negative  gravity  points  downward 

y — cp/c,  - ratio  of  specific  heats  of  the  fluid 

K » thermal  conductivity  coefficient  (erg/cm-sec-°K) 

h,  — enthalpy  of  the  ith  species  of  the  fluid  (ergs/molecule) 

V’,  — diffusion  velocity  of  the  ith  species  (cm/sec) 
v — viscosity  coefficient  (gm/cm/sec) 

P — pressure  (dynes/cm2) 
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V — volume  of  the  container  (cm3) 

Du  - binary  diffusion  coefficient  for  species  "i"  through  species  "j"  (cm2/sec) 

Si  — diffusion  source  term  for  the  ith  species  (cm-1) 

As  this  is  a two-dimensional  calculation,  both  the  (/)  and  (z)  components  of  V and  V'  are  car- 
ried. In  addition,  we  require  two  more  equations  to  close  this  set.  We  have  used 

P - nkT  (6) 

for  the  relation  between  pressure,  density  and  temperature,  and 


V,J 


to  relate  the  diffusion  velocities  to  the  other  physical  parameters.  The  quantities  S,  and  A„  are 
defined  and  discussed  in  detail  in  Appendix  1.  The  quantities  IS,}  are  forcing  terms  for  molecu- 
lar diffusion  which  include  density,  pressure,  and  temperature  gradients  as  well  as  differential 
body  forces  on  the  various  species. 


This  2D  model  has  a rather  general  initialization  and  diagnostics  capability  in  order  to  con- 
sider a variety  of  problems.  The  numerical  model  is  designed  to  operate  efficiently  on  a vector 
computer,  thus  demonstrating  the  value  of  vector  computation  even  in  complex  problems  such 
as  reactive  flow. 


This  paper  discusses  some  of  the  details  of  the  model,  that  is  the  actual  form  of  Eqs.(l-7) 
used,  the  simplifying  assumptions  made,  and  their  justifications.  In  Section  II  areas  of 
potential  application  and  the  corresponding  numerical  difficulties  which  they  present  are  dis- 
cussed. Section  HI  describes  the  numerical  model  in  greater  detail.  The  grid  system,  spatial 
differentiation,  timestepping,  diffusive  transport,  and  chemical  kinetics  algorithms  are  discussed. 
Section  (IV)  discusses  calculations  which  have  been  performed  to  benchmark  the  model  in 
several  difficult  but  well  understood  situations.  These  tests  include  comparisons  with  a one- 
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dimensional,  fully  compressible  Lagrangian  code  and,  for  simple  cases,  comparison  with  a one- 
point  analytic  calculation.  Agreement  is  quite  good,  and  we  have  been  able  to  define  the  limits 
of  accuracy  of  the  approximations  and  numerical  algorithms  which  are  used. 

II.  Problems  to  be  Solved 

There  are  several  types  of  cylindrically  symmetric  problems  to  which  this  model  can  be 
applied.  One  example  is  a laminar  gas  jet,  such  as  a bunsen  burner  or  a flat  plate  burner,  where 
the  fuel  and  oxidizer  may  mix  either  prior  to  or  after  entering  the  combustion  chamber.  In 
either  case  there  is  no  reason  to  expect  macroscopic  asymmetry  in  the  azimuthal  direction. 
However,  species  densities  and  temperatures  very  dramatically  both  above  the  ignition  point 
and  radially  outwards  from  the  jet  center,  even  in  steady  state.  The  two  problems  differ 
significantly  in  the  initial  conditions  required  and  the  way  in  which  the  geometric  boundary  con- 
ditions are  applied.  One  is  a diffusion  flame  and  the  other  a premixed  flame.  In  the  case  of  the 
premixed  system,  the  processes  of  interest  will  be  the  jet  and  flame  speed  and  the  detailed 
chemical  kinetics.  For  the  diffusion  flame,  however,  the  physics  is  made  more  complicated  by 
the  fact  that  a finite  time  is  required  for  the  oxidizer  to  mix  with  the  fuel.  Computational 
models  which  are  to  be  applied  to  these  typical  practical  situations  invariably  require  a non- 
uniform  grid  to  resolve  the  combustion  regions  finely  while  limiting  the  overall  number  of  grid 
points  at  which  computations  must  be  made.  The  problem  here  is  one  of  widely  disparate  physi- 
cally important  space  scales. 

Buoyancy  is  a major  effect  in  most  diffusion  flames  and  even  in  premixed  systems  it  can 
be  responsible  for  much  of  the  flame  dynamics.  The  computational  problem  introduced  here  is 
one  of  widely  disparate  fluid  dynamic  timescales.  The  slow  flow  algorithm  used  for  the  tem- 
poral integration  of  the  fluid  dynamic  equations  is  a way  to  filter  sound  waves  out  of  the  equa- 
tions so  that  timesteps  much  longer  than  the  Courant  step  8 1 — 8x/Cs  can  be  taken  without 
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numerical  instability.  The  standard  approach  of  formulating  an  implicit  pressure  equation 
requires  at  least  formal  substitution  of  one  equation,  in  finite  difference  form,  into  another  and 
usually  some  form  of  additional  numerical  smoothing.  The  resulting  algorithms  do  not  resolve 
steep  gradients  well  during  convection.  Flux-Corrected  Transport  (FCT)  is  a far  superior  con- 
vection algorithm  but  is  intrinsically  nonlinear  and  hence  does  not  lend  itself  to  implicit  formu- 
lations. Therefore  the  slow  flow  algorithm,  which  is  asymptotic  rather  than  implicit  in  concept, 
allows  the  unlimited  use  of  FCT  to  convect  steep  species  and  temperature  gradients.  By  filter- 
ing out  the  fast  sound  waves,  the  physically  important  buoyancy  and  molecular  diffusion  times- 
cales can  be  integrated  inexpensively  using  only  a few  timesteps  rather  than  hundreds  of 
thousands.  ' * ' 

Another  class  of  important  problems  concerns  turbulent  .nixing  with  reactions.  Large- 
scale,  quasi-s'ationary  eddy  structures  are  observed  in  most  shear  flow  situations  and  these 
structures  dearly  play  a vital  role  in  the  turbulent  mixing  process.  While  the  full  details  of  fully 
developed  turbulent  flow  is  beyond  numerical  simulation,  these  eddy  structures,  which  are 
often  two-dimensional,  can  be  simulated  accurately.  Therefore  FLAME  is  designed  to  allow 
calculation  of  these  structures  with  and  without  self-consistent  energy  release  from  the  chemical 
kinetics.  The  importance  of  convective  acceleration  in  this  kind  of  problem  necessitates  carry- 
ing the  slow  flow  expansion  to  first  order  in  the  perturbed  pressures,  but  this  is  relatively  easy 
to  do.  The  goal  is  to  have  a fluid  dynamics  model  which  can  resolve  these  structures  realisti- 
cally as  a basis  for  testing  macroscopic  turbulence  closure  models. 

Another  area  of  applicaton  for  FLAME  is  the  phenomenological  modelling  of  fire  spread 
in  enclosed  spaces.  For  this  class  of  problem,  the  extreme  care  taken  to  provide  an  accurate 
fluid  algorithm  with  good  timestep  properties  is  almost  overkill.  Nevertheless,  since  FLAME  is 
designed  for  detailed  modelling,  the  extra  accuracy  and  reliability  of  the  fluid  dynamics  should 
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pay  dividends  as  a basis  for  phenomenological  and  empirical  modelling  as  well.  To  this  end, 
methods  for  representing  inflow-outflow  conditions  and  interior  structures  and  obstacles  are 
being  studied.  Before  such  a model  becomes  complete,  however,  extensive  additional  numeri- 
cal work  will  have  to  be  undertaken  to  build  a multidimensional  radiation  transport  model  capa- 
ble of  dealing  equally  accurately  with  volumes  which  are  optically  thick  as  well  as  optically  thin. 

III.  Numerical  Model 

There  are  three  aspects  to  the  integration  of  the  equations  which  describe  the  system:  1) 
spatial  integration  and  differentiation;  2)  time-stepping  and  integration  of  the  energy  (pressure) 
equation  using  the  assumption  of  asymptotic  relaxation,  i.e.  "slow-flow";  3)  physical  transport 
and  chemical  kinetics. 

III.l.  Grid  System  — Spatial  Integration  and  Differentiation 

The  cylindrical  finite  difference  grid  used  in  the  numerical  model  is  shown  in  Fig.  (1). 
Spatial  integration  and  differentiation  is  performed  on  this  staggered  grid.  The  representation  of 
the  fluid  dynamics  is  inherently  Eulerian,  although  rezoning  and  variable  grid  spacing  are 
allowed  at  each  time  step.  This  allows  computation  of  quantities  such  as  V • V to  fall  automati- 
cally on  the  grid  points  where  they  must  be  used.  For  example,  for  the  continuity  equation, 

^2.  = -p^7-  V,  whereas  p is  defined  on  the  grid,  V must  be  defined  at  interstitial  points  in 
dl 

order  for  V V to  appear  on  the  same  grid  as  p.  In  solving  for  V,  the  reverse  calculation  of  the 
velocity  from  and  ^ and  <J>,  the  stream  function  and  velocity  potential,  is  similarly  straightfor- 
ward if  <//  and  <t>  are  defined  at  the  proper  places.  When  mass  flux,  that  is  momentum,  is  a 
necessary  quantity  such  as  in  calculating  the  diffusion  velocities,  the  only  quantities  which  then 
need  to  be  averaged  are  the  mass  densities  at  the  half-cell  positions.  Thus  this  staggered  grid 
structure  insures  that  effectively  second-order  accuracy  is  maintained  in  spatial  integration,  at 
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least  for  p uniform  or  slowly  varying  grid  spacing.  The  grid  and  numerical  model  in  general  can 


be  applied  to  Cartesian  and  spherical  as  well  as  cylindrical  coordinates.  The  model  is  set  up  to 


make  the  change  relatively  easy,  but  the  geometrical  constants  used  in  spatial  differentiation  as 


well  as  the  axis  boundary  condition  (AT?  = 1)  would  have  to  be  changed 


The  quantities  defined  at  grid  points  are  D,  £ p,  (p,|,  T.  p and  d>.  The  stream  function 


ip  and  velocity  V are  carried  at  the  interstitial  points  shown  in  Fig.  (1).  The  grid  is  staggered  in 


this  manner  to  1)  facilitate  calculation  of  V from  <l>  and  ip  (and  the  reverse  of  course),  2)  to 


make  the  model  conservative  without  resorting  to  exotic  boundary  conditions  and,  3)  to  sim 


plify  the  implementation  of  boundary  conditions.  In  order  to  insure  conservation  of  mass,  for 


example,  in  a dosed  system,  the  flux  of  mass  across  a boundary  must  be  zero.  Thus  we  choose 


0 and  ip  = 0 at  the  outside  cylindrical  wall,  as  is  shown  in  Fig.  (1) 


The  self-consistent  nature  of  this  grid  can  best  be  seen  by  use  of  an  example.  Consider 


the  V ■ V term  in  Eq.  (1)  as  a source  term  in  the  Poisson  equation 


The  term  V V becomes 


The  term  becomes  (for  equal  spacing  in  R and  Z) 


tion  specified  at  the  same  physical  place.  A similar  relationship  exists  for  V and  i<i  (or  V x V 
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and  V X V X but  with  values  for  $ on  a different  staggered  grid  than  that  used  for  <J>. 
Thus  the  four  velocity  variables  ♦,  i|r,  and  Kz  must  be  on  different  grids  in  order  to  be 
self-consistent  and  conservative  without  i sorting  to  artificial  boundary  conditions  and/or  excess 
averaging. 


Equations  (2)  and  (3),  the  curl  and  divergence  of  velocity  respectively,  require  the  calcu- 
lation of  velocity  itself  as  a source  term.  In  addition,  velocity  is  a useful  diagnostic  tool  and  is 
used  for  other  calculations  as  well.  For  simple  theories,  or  low  order  approximations,  the  vorti- 
city  and  divergence  alone  could  be  used.  Such  is  not  the  case  for  our  model  and  hence  we 
need  to  be  able  to  calculate  the  velocity  from  VxV  and  V • V.  This  is  easily  done  by  defining 
a stream  function  and  velocity  potential  4>  by 

V - V<t>  + Vx«jT. 

This  then  leads  to  two  equations  for  4>  and  i /». 


V2d>  - V • V - D 


(10a) 


(10b) 


(11) 


and 

VxVx^  - VxV  = f . 

These  are  elliptic  equations  of  the  general  form 

A + B -fj  + C -r^r-  + />■—  + £“■  + G|* 
dx2  By2  BxBy  Bx  By 

with  either  Neumann  or  Dirichlet  boundary  conditions,  both  of  which  arise.  Part  of  our  origi- 
nal choice  for  the  staggered  grid  was  based  on  being  able  to  set  these  boundary  conditions 
easily.  Thus  Eq.  (10)  can  be  cast  in  the  form 

V2d>-‘F,  (12a) 


and 


VxVx*-f2  (12b) 

which  are  elliptic  equations  of  the  general  form  of  Eq.  (11).  The  substitution  method  of 

Madala2  solves  these  equations  more  quickly  than  the  usual  relaxation  methods  and  with  more 
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accuracy  than  iterative  schemes.  Most  important,  however,  is  that  it  allows  the  coefficients  A, 
B,  C,  D,  E and  F to  be  functions  of  both  position  and  time  which  allows  for  stretched  grids  and 
simole  rezoning. 

III. 2.  Time  Stepping  and  the  Asymptotic  Approximation 

We  have  incorporated  time-splitting  in  order  to  make  the  model  flexible  and  to  simplify 
adding  modules  to  describe  specific  physical  effects.  Except  for  the  driving  program,  which 
invokes  the  hydrodynamics,  molecular  diffusion  and  thermal  conductivity,  the  various  modules 
which  incorporate  the  physical  processes  are  exercised  separately  and  interact  as  source  terms  at 
each  time  step.  This  splitting  is  a standard  technique  but  the  inherent  assumptions  must  always 
be  tested  when  implementing  a new  numerical  model. 


In  FLAME  all  time-stepping  is  done  explicitly  in  a time-centered  fashion.  Doing  the  cal- 
culation explicitly  allows  more  flexibility  in  the  model  development.  Time-centering  ensures 
reversibility  in  the  physics.  By  reversing  the  sign  of  the  time  step,  8r,  we  can  make  the  reversi- 
ble physics  run  backwards.  This  is  a real  property  in  non-dissipative  physical  systems  which 
ought  to  be  mirrored  as  closely  as  possible  in  a numerical  model.  When  diffusion  is  included, 
such  reversibility  ceases  to  be  strictly  valid  of  course. 


Although  explicit  methods  have  time  step  limitations  to  avoid  instability,  our  asymptotic 
pressure  relaxation1  algorithm  is  an  exception  to  this  general  rule.  To  illustrate  this  point,  we 
note  that  Equations  (1),  (2),  and  (5)  are  solved  as  is.  However,  Equations  (3)  and  (4)  are 
combined  to  form  the  following  equation 


(13) 
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where 


chem 


Such  a formulation  embodies  the  assumption  that  waves  travelling  at  the  sound  speed  do  not 


influence  the  overall  hydrodynamics  motion.  This  is  the  reason  for  the  name  "Slow-Flow1 


This  formulation  excludes  shock  waves  from  the  set  of  problems  which  we  can  consider 


The  pressure  can  still  vary  overall,  but  now  only  on  the  timescale  of  the  chemical  energy 


release.  This  is  the  time-stepping  approach  used  in  the  solution  of  the  energy  equation.  For 


Eqn.  (1),  (2)  and  (5)  which  are  solved  explicitly,  we  use  flux-corrected  transport 


method  is  both  relatively  non-diflusive,  and  very  fast.  The  worst  problem  which  arises  using 


this  method  is  the  non-linear  clipping  phenomenon,  which  will  be  discussed  in  detail  in  the 


next  section  under  results.  Generally  clipping  occurs  when  there  is  a sharp  peak  in  the  mass 


density  (p)  or  vorticity  (£).  The  anti-diffusive  flux  terms  (which  allow  very  steep  gradients  to 


propagate  conservatively)  tend  to  flatten  such  peaks.  As  will  be  seen,  this  introduces  little  error 


even  in  the  worst  case,  e.g.,  a Heaviside  function  of  density  or  vorticity  perturbation,  and  the 


error  is  virtually  non-existent  when  small  gradients  are  present 


The  numerical  timestep  to  which  we  are  limited  is  the  shortest  of  those  required  to 


integrate  (stably)  molecular  diffusion,  thermal  conductivity  and  fluid  flow.  Chemistry  computa 


tion  is  allowed  to  subcycle  itself,  so  this  will  not  limit  the  overall  timestep  unless  so  much 


energy  is  released  that  it  affects  the  fluid  flow  and  diffusion  timestep  limits.  Typical  time  con 


slants  for  molecular  diffusion  can  be  found  from  the  diffusion  equation  cast  into  the  form 


V (DVp) 
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with 


D = 


1QI7T3/4 

n 


(17; 


Thus,  for  example,  for  temperatures  in  the  range  of  500-1000°K  and  a grid  spacing  of  0.25  cm, 
the  diffusion  limit  timesteps  will  be  about  0.1-0.01  seconds.  These  limits  vary  with  time  and 
must,  of  course,  be  checked  at  each  timestep. 


The  equations  are  integrated  using  a time  centered  scheme.  The  order  of  the  integrations 
is  as  follows: 


1)  Estimate  values  for  P,  p,  {p,},  £,  D,  T,  i /»,  at  d>  at  the  new  time. 

2)  Find  source  terms  for  Eqn.  (1,  2,  5 and  10)  based  on  time  centered  values 
(/„  + 1/2  A/). 

3)  Integrate  Eq.  (1,  2,  5 and  10)  from  (r„)  to  (/„  + A t)  using  the  source  terms 
defined  at  the  time  centered  position.  This  is  a fluid  flow  calculation  only. 


4)  Repeat  steps  (2)  and  (3)  until  convergence  is  reached. 

5)  Integrate  the  kinetics  portion  of  Eqn.  (5,  6,  13,  14  and  15)  using  the  Asymptotic 
Chemical  Kinetics  Scheme  of  Young  and  Boris.6  The  subcycling  done  at  this 
point  effectively  integrates  the  equation  of  siate,  and  the  energy  release  is  accu- 
mulated as  a source  term  for  step  (1). 

The  following  illustration  shows  where  each  of  these  time-step  points  is,  in  relation  to  steps  (1-5) 

tn  r„  + A l 

1)  • =>•  estimate  values  at  (/„  + A t) 

to  to  + At 


2)  * > X < • find  sources  at  ( t . + A 1 
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I 

i 

l 


3) 

4) 

5) 


t„  + 1/2  Af 


integrate  from  (r„)  to  (r„  + At) 
repeat  steps  (2  and  3) 

integrate  chemical  kinetics 


Since  each  step  of  the  integration  is  of  at  least  second  order  accuracy,  the  overall  scheme  will 
also  be  second  under  accurate.  In  addition,  using  the  time  centered  source  terms  in  steps  (2 
and  3)  insures  reversibility,  at  least  for  non-diffusive  physical  systems,  once  the  iterative 
scheme  has  converged. 


1II.3.  Physical  Transport  and  Chemical  Kinetics 

There  are  two  pieces  of  physics  built  into  the  model  which  are  crucial  to  model  any 
combustion  problem  accurately.  Molecular  diffusion  and  thermal  conductivity  have  been  built 
into  the  basic  model.  Both  are  essentially  source  terms  in  Eq.  (13).  At  present,  we  use  the 
conceptually  simple  results  of  Sears7  for  KT,  namely 

K'-myPxy/r.-  (18) 

where  A is  the  average  molecular  mass  of  the  gas  and  T„  — 293°K.  For  molecular  diffusion, 
we  calculate  the  actual  diffusion  velocities  for  each  species  from  the  fully  coupled  diffusion 


equations.  The  method  is  discussed  in  detail  in  Appendix  I 


The  chemistry  portion  of  our  model  is  a general  five-species  model.  Such  a configuration 


was  chosen  so  that  we  could  model  such  systems  as  C0-02  and  Oj-02  in  detail,  or  H2-02  and 


and  CH4-02  using  global  schemes.  Appendix  II  shows  a simple  model  scheme  and  also  the 


one  used  for  C0-02  bum.  When  chemical  kinetics  is  present,  Eq.  (3)  must  be  solved  for  crea 


tion  and  consumption  of  individual  species.  As  was  done  for  the  other  physical  effects,  energy 


release  or  absorption  is  a source  term  in  the  energy  equation,  Eq.  (13) 





mmtm 


m 
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As  in  our  choice  of  grid  and  hydrodynamic  algorithms,  we  must  strike  a balance  between 
physical  complexity  and  model  flexibility.  More  species  can  be  added,  of  course,  but  having 
five  species  allows  a large  number  of  important  calculations  to  be  performed  without  detracting 
from  the  main  purpose  of  FLAME  which  is  to  perform  accurate  multidimensional  calculations. 
With  five  species  we  can  model  the  C0-02  and  the  03-02  systems  in  detail.  The  H2-02  sys- 
tem and  more  complex  hydrocarbon  systems  such  as  methane  and  benzene  can  be  modelled 
using  global  schemes.  For  the  simple,  detailed  models  we  can  even  add  diluents  to  test  for 
flammability  and  flame  extinction.  For  the  global  models  of  reactive  systems  the  emphasis 
shifts  to  studying  detailed  fluid  dynamic  effects.  The  flow  of  trace  species  in  nitrogen  pressuri- 
zation experiments  can  be  followed  and  more  realistic  reaction  schemes  than  A + B — C can 
be  used  for  turbulent  mixing  calculations  and  calibration  of  turbulence  closure  models. 


IV.  Tests  of  the  Model 


In  order  to  be  able  to  trust  the  predictive  capability  of  the  model,  as  a minimum  we  must 


insure  that  it  behaves  properly  in  cases  where  we  know  the  answer.  The  two  aspects  to  check- 
ing this  predictive  capability  involve  benchmarking  it  against  well  documented  and  thoroughly 
tested  numerical  models  and  also  against  experiments  for  which  we  can  reasonably  expect  good 
agreement  and  which  have  been  well  diagnosed. 


We  report  here  on  the  former  set  of  tests,  namely  comparison  with  other  models  and  cal- 
culations. In  each  case,  where  possible,  we  have  compared  this  model  with  a one-dimensional 
fully  implicit  Lagrangian  code  and  also  a simple  "analytic"  model.  The  tests  we  applied  were 
both  passive,  that  is  no  reactions  occurred,  and  active  and  hence  included  reactions.  The 
former  tests  check  the  internal  consistency  of  the  hydrodynamics  whereas  the  latter  check  the 
validity  of  the  coupling  of  chemistry  to  the  hydrodynamics.  In  both  sets  of  tests,  we  used  the 
Heaviside  temperature  perturbation8  shown  in  Fig.  (2)  and  given  by 
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Po  - Ap,  r < )?0| 
p„  ,>*„(  *■ ' J 67cm 

The  first  set  of  tests  were  passive,  that  is  no  reactions  were  present.  This  would  be  the 
situation  if  only  one  species  were  present,  or  if  the  two  initial  species  were  oxygen  and  nitro- 
gen. All  combinations  of  molecular  diffusion  and  conductivity  were  tried.  Figure  (3)  shows 
the  comparison  between  FLAME  (•)  and  the  Lagrangian  one  dimensional  model  (O)  after 
three  seconds  with  both  mass  diffusion  and  thermal  conductivity  included.  By  plotting  the  tem- 
perature as  a function  of  radius  in  these  figures  the  2D  data  are  collapsed  into  a simple  ID  plot. 
The  spread  or  scatter  of  the  various  FLAME  temperature  values  near  a given  radius  gives  a 
good  measure  of  the  numerical  truncation  error  because  this  scatter  arises  from  performing  a 
spherical  calculation  on  a cylindrical  r-z  mesh.  As  can  be  seen,  the  agreement  between  the  two 
is  quite  good  even  for  relatively  long  times.  Also,  there  is  very  little  asymmetry  in  the  2D  cal- 
culation. Near  the  peak,  however,  there  is  some  clipping  due  to  the  use  of  the  FCT  routines  on 
the  rather  coarse  grid  used  for  these  tests  (20x20).  To  assure  ourselves  that  indeed  the  two 
codes  would  agree  qualitatively  we  did  one  test  on  a (40x40)  grid,  shown  in  Fig.  (4).  Since  we 
see  errors  of  less  than  5%  in  both  cases  the  general  problem  of  convective  transport  will  be  well 
handled,  even  on  the  rather  coarse  grid.  When  we  turn  to  the  problem  of  turbulence,  however, 
we  will  need  the  resolution  at  particular  grid  points  and  then  we  can  either  use  fine  gridding,  or 
variable  grid  spacing  to  obtain  high  resolution  at  the  necessary  points.  The  latter  technique, 
called  variable  rezoning,  is  presently  being  worked  on. 

One  necessary  test  was  a check  of  molecular  diffusion  when  only  one  species  was  present, 
that  is 

P “ Pa- 

ox  when  equal  parts  of  two  species  are  present,  which  have  the  same  mass,  and  are  in  all  other 
aspects  identical.  In  such  a situation,  there  should  be  no  net  diffusion.  This  is  the  same  as  say- 
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j 

j 

I 

I 
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ing  that  there  is  no  diffusion  of  the  total  mass  density  due  to  molecular  diffusion.  This  was 
verified  for  both  the  one  and  two-dimensional  models  and  is  attributable  to  the  manner  in 
which  we  calculate  the  diffusion  velocities. 

These  tests  have  validated  the  internal  consistency  of  the  hydrodynamic  parts  of  the 
model.  Next  we  tested  the  basic  assumption  of  time-step  splitting.  The  first  check  was  of  the 
chemistry  in  a configuration  where  hydrodynamics  could  not  play  a part.  Such  is  the  case  when 
uniform  ignition  is  used,  namely, 

P “ Pa  + Pb • Pa  ” Pb 
and 

p - p„  - Ap 

for  all  radii.  For  these  conditions  there  will  be  no  transport.  All  three  models  (including  the 
one  point  analytic  model)  agreed  to  within  one  degree  Celsius,  which  is  adequate.  For  asymp- 
totic calculations  of  the  sort  in  which  we  are  interested,  which  last  for  tens  of  seconds,  shorten- 
ing the  timestep  to  the  chemistry  time  scale  (as  was  necessary  in  the  previous  calculation)  is 
not  a satisfactory  solution.  Instead,  we  chose  to  subcycle  the  chemical  kinetics  and  couple  this 
to  the  hydrodynamics  by  accumulating  the  source  terms  for  all  subcycles  and  updating  the  equa- 
tion of  state  at  each  subcycle.  An  independent  equation  for  the  temperature  based  on  energy 
release  could  have  been  solved,  but  we  have  found  stability  problems  in  this  approach.  The 
technique  of  subcycling  allows  us  to  use  a hydrodynamics  timestep  of  one  to  two  orders  of  mag- 
nitude larger  than  that  required  for  the  chemical  kinetics. 

The  validity  of  this  approach  is  demonstrated  by  the  calculation  which  is  shown  in  Fig. 
(5).  In  (Sa)  the  bubble  radius  is  plotted  and  in  (Sb)  the  peak  temperature  is  shown,  both  as  a 
function  of  time.  In  both  cases,  the  points  (O)  represent  the  (2D)  calculations,  and  the  points 
(•)  represent  the  equivalent  (ID)  calculation.  As  can  be  seen,  the  peak  temperatures,  ignition 


15 
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times,  and  bubble  radius  agree  satisfactorily.  For  these  computations,  only  chemistry  and  the 


fluid  flow  are  allowed  to  play  a part.  Both  molecular  diffusion  and  thermal  conductivity  are 


absent 


To  complete  this  series  of  tests,  which  checks  the  symmetry  of  the  two-dimensional 


model  as  well  as  comparing  it  with  a well  documented  one-dimensional  model  and  a single 


’analytic”  one  point  model,  we  ran  a test  with  molecular  diffusion,  thermal  conductivity,  client 


istry  and  hydrodynamics,  all  interacting.  Figure  (6)  shows  a plot  of  the  temperature  profile  for 


FLAME  (•)  and  the  Lagrangian  one-dimensional  model  (O),  in  this  configuration.  As  can  be 


seen,  differences  between  the  three  models  are  quite  small.  Temperature  has  been  used  as  a 


diagnostic  for  two  reasons:  first,  it  is  sensitive  to  errors  in  both  the  hydrodynamics  and  chemi 


cal  kinetics;  second,  it  is  a straightforward  physical  quantity  and  is  easily  visualized  and  meas 


ured.  Asymmetries  which  appear  can  then  be  discussed  in  terms  of  too  great  (or  too  small)  a 


flow  or  too  much  energy  release.  If  the  two  dimensional  hydrodynamics  were  done  incorrectly, 


it  would  most  likely  show  up  as  an  asymmetry,  and  if  either  the  hydrodynamics  or  the  chemis 


try  were  incorrect,  a difference  between  the  temperature  profiles  of  the  two  models  would  be 


apparent 


As  can  be  seen,  the  FLAME  implementation  of  the  reactive  flow  model  in  two  dimen 


sions  agrees  well  with  both  one  point  analytic  calculatons  as  well  as  with  a well  tested  one 


dimensional  model.  Furthermore,  it  is  conservative  both  for  mass  and  energy.  The  early  prob- 


lems to  be  modeled  should  be  those  for  which  well  diagnosed  experiments  exist.  Eventually 


however,  one  will  want  to  attempt  predictions  for  systems  which  are  in  the  design  phase  or  for 


which  experimental  data  are  difficult  or  impossible  to  obtain 
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The  potential  applications  of  this  model  are  numerous  and  expand  still  further  when  cer 


tain  modifications  and  extensions  of  the  model  are  contemplated.  Using  the  code  in  its  current 


form,  a number  of  gas  jet  and  burner  problems  are  being  considered  which  will  contribute  to 


the  interpretation  of  chemical  kinetics  investigations  ongoing  in  a number  of  laboratories.  In 


these  cases  rather  detailed  kinetics  packages  will  have  to  be  incorporated  but  laminar  diffusion 


and  premixed  flames  will  be  used  so  that  turbulence  will  not  cloud  the  reaction  kinetics  issues 


The  basic  FLAME  code  will  be  equally  useful  in  treating  heterogeneous  combustion  and  tur 


bulent  mixing  because  the  basic  hydrodynamics  is  about  as  accurate  and  non-diffusive  as  is  pos 


sible  in  an  Eulerian  representation.  By  simplifying  the  chemistry  and  concentrating  on  the  flow 


complications,  extensive  turbulent  fields  can  be  simulated  with  and  without  subgrid  turbulence 


closure  models.  In  special  cases  FLAME  can  be  used  to  evaluate,  via  simulation,  the  probabil 


ity  density  functions  which  are  often  invoked  in  turbulence  modelling,  but  are  not  known  accu 
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Appendix  I 


Multi-species  diffusion  velocities  can  be  related  to  those  due  to  binary  diffusion  by  consid 


ering  the  general  diffusion  equation1 


V(«,/n)  - £ 


where  N is  the  number  of  species,  pl7  is  the  reduced  mass,  //  represents  the  body  force  on  the 


/ th  species  and  vti  is  the  binary  viscosity  coefficient.  If  we  use  the  usual  definition  for  the 


binary  diffusion  coefficient,  namely 


then  Eq.  (1.1)  becomes 


We  can  cast  this  into  an  implicit  form  by  defining  as  a source  term,  all  those  elements  of 


Eq.  (1.3)  which  are  not  explicitly  velocity  dependent,  namely 


S,  = V(njn) 
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The  third  source  term  in  Eq.  (1.4)  is  diffusion  due  to  thermal  gradients.  It  is  included  because 


there  are  situations  in  which  it  can  be  important,  especially  when  there  are  heavy  and  light 


species  mixing,  e.g.,  hydrogen  and  oxygen.  With  this  definition  for  source  terms,  Eq.  (1.3)  can 


be  written  in  the  form 


Since  -0  and  on  the  right  hand  side  of  Eq.  (1.5)  there  are  only  N-  1 independent  terms 


Eq.  (1.5)  is  singular.  Therefore,  we  require  one  more  equation  to  close  this  set.  For  this,  we 


use  the  physical  condition 


which  implies  that  the  total  mass  continuity  equation  has  no  net  diffusion  term 


In  order  to  make  this  a more  tractable  problem,  we  assume  that  the  diffusion  velocities 


separate  into  directions  associated  with  the  vector  gradients.  Then  the  equations  to  be  solved 


become  two  dimensional  matrix  equations  (in  species),  rather  than  a full-blown  tensor  equa- 


with  the  mutual  term  dropped  explicitly.  This  term  allows  us  to  define  a diffusion  coefficient 


for  each  species  relative  to  the  rest  of  the  fluid,  in  the  spirit  of  the  binary  diffusion  coefficient 


Define 


as  an  average  diffusion  coefficient.  It  has  the  following  properties: 
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1.  In  the  case  of  two  species,  it  reduces  to  the  binary  diffusion  coefficient 


2.  It  has  the  same  symmetry  properties  as  the  reduced  mass  (/i,;) 


3.  It  will  be  the  same  for  two  species  which  behave  identically 


4.  Reduces  to  zero  for  a single  species. 


S.  Allows  us  to  make  a perturbation  expansion  for  the  diffusion  velocities  in  terms  of 


the  diffusion  coefficient  and  the  source  term,  S, 


Defining  a perturbation  series  for  V\  we  have 


and  substituting  into  Eq.  (1.8)  yields 


A judicious  choice  for  Vf 


will  allow  us  to  make  the  corresponding  expansion  for 


By  substitution  we  find 


etc.  where 


Combining  terms  yields,  given  the  above  definition  for  A. 


where  the  vector  S is  now  over  species.  The  problem  we  a.c  left  with  is  the  convergence  of 
this  series.  The  convergence  criterion  is  that 


Numerically,  however,  there  may  be  a problem  with  differences  of  large  numbers  produc 


ing  round-off  errors.  Addition  of  a constant  (C)  to  each  row  of  matrix  A does  not  change  the 


value  of  8 V„  since  by  application  of  Eq.  (1.16),  the  transformation 


A )k  “ Alk  — C,  (1.17) 

leaves  6 V,,  68  V,  etc.,  unchanged.  Since  the  choice  of  (Cj  is  arbitrary,  we  chose  it  so  that  the 


mean  value  of  each  row  of  A vanishes 


and  thus  the  "difference-in-large-numbers"  problem  is  solved  prior  to  finding  the  8 l^s  etc. 
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We  have  implemented  the  above  diffusion  model  in  a subrountine  called  DFLUX,  which 


has  been  optimized  for  a vector  computer.  In  the  two  dimensional  model,  we  find  the  diffusion 


fluxes  for  all  species  for  each  layer  in  the  radial  direction.  Thus  the  natural  vector  length  will 


be  NR  x A fSP  in  length  where  NR  represents  the  number  of  grid  points  in  a horizonal  layer 


and  NSP  is  the  number  of  species.  Figure  (7)  shows  the  time  required  to  calculate  the  {/>,  V,} 


species  at  each  grid  point.  As  can  be  seen,  the  computation  time  required  for  each  grid 


increases  much  more  slowly  than  the  number  of  grid  points.  A more  dramatic  demonstration 


of  this  is  shown  in  Figure  (8)  where  the  time  required  for  each  grid  point  is  plotted  as  a func 


tion  of  the  layer  size,  with  the  total  grid  being  held  constant  at  400  points.  The  times  are  nor' 


malized  to  the  time  required  for  a single  point  (+)  and  for  a vector  of  length  ten  (•).  The  two 


curves  approach,  asymptotically,  the  speed  limitation  of  the  computer  hardware. 


For  more  than  about  ten  grid  points  this  method  is  faster  than  calculation  of  the  usual 


diffusion  source  term  V • (D  7 p)  usually  used  as  a source  term  in  the  momentum  equation 


and  in  addition  is  a more  satisfactory  approach 


Appendix  II 

Both  the  model  chemistry  and  the  carbon  monoxide  schemes  are  five  species  models.  We 


have  used  the  model  chemistry  for  testing,  however,  since  we  know  its  time  history  on  a 


theoretical  basis. 


For  the  model  chemistry,  species  four  and  five  are  used  simply  as  diluent.  This  allows  a 


more  satisfactory  treatment  of  the  scheme  A + B — ► C due  to  momentum  and  energy  conser 


vation  considerations,  but  is  solved  just  as  quickly.  All  rates  are  in  a simplified  Arrhenius  form 


R,  - C,  exp  (-  EJT). 

The  values  for  the  coefficient  C,  and  E,  are  given  in  Table  (I) 

Table  I 


The  scheme  is 


This  scheme  thus  allows  for  both  two  and  three  body  interactions. 


c 

E 

4.8  x 10",J 

6.0  x 10J 

5.0  x 10“w 

6.0  x 10J 

9.6  x 10+4 

7.0  x 10J 

1.0  x 10'14 

7.0  x 10* 

°J. 

7.0  x 10J 
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This  carbon  monoxide  scheme  is  somewhat  more  complex.  We  solve  the  following  set  of 
equations  which  also  correspond  to  the  source  terms  in  Eq.  (S). 

-^ICO]  - -C2/[C0][0][M]  --  C3/lCO][Oj  01.4) 

+ CjrICOjHM]  + C3*[C02][0] 

OJ  - + C , y[0] 2[M]  - C3/IC0][02]  (11.5) 

-C,*[02][M]  + C3*[C02][0] 

■j[0) 2Ct/[0] 2(M]  - C2/[C0110UM]  + C3/[C0][02]  01.6) 

+ 2C,„[02][M]  + C2*[C02][M)  - C3*[C02][0] 

-^[C02]  - + C2/[COHOHM]  + C3/[COU02l  (11.7) 

- C2*[C02HM]  - C3r[C02][01. 

We  then  have  the  following  correspondence  with  the  {/»>}, 


«,  - [CO] 

n2  - lo2] 

«3  - [o] 

*4  - [C02] 


dn  5 

and  n5  — a diluent  whose  total  density  is  fixed,  that  is  = 0.  These  equations  correspond 


to  the  following  reactions 


Forward  (f) 

Reverse  (r) 

(1) 

O + O + M — Oj  + M 

02+M  — O + O + M 

(2) 

CO  + 0 + M — C02  + M 

C02  + M — CO  + 0 + M 

(3) 

CO  + 02  — C02  + 0 

C02  + 0 - CO  + 02 . 

The  general  form  of  the  rate  coefficients  is 

C - C^Tae~Err 

and  the  coefficients  are  given  in  Table  (II). 
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OUTSIDE  WALL  (NR  + I ) 


CYLINDER  AXIS 


Schematic  diagram  of  the  grid  system  used  in  the  two-dimensional  flame  model.  The  figure 
shows  a staggered  grid  and  indicates  where  each  variable  is  defined. 
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INITIAL  TEMPERATURE  PROFILE 


0.0  RADIUS  (CM)  | 

Fig.  2.  The  initial  temperature  distribution  used  in  the  calibra 
lion  calculations  or  the  two-dimensional  model.  The  tempera 
lure  is  plotted  as  a function  of  radius. 


0.0  RADIUS  (CM)  ► 

Fig.  3.  Scatter  plot  of  temperature  verus  radius  when  both 
molecular  diffusion  and  thermal  conductivity  are  included.  The 
scatter  arises  from  doing  a spherically  symmetric  calculation  on  a 
cylindrical  grid.  The  initial  condition  for  this  calculation  is 
shown  in  Fig.  (2).  A scatter  plot  of  this  type  reveals  any  grid  in- 
duced asymmetries  in  the  2D  numerical  model.  In  addition,  this 
representation  is  a convenient  vehicle  for  comparison  with  the 
ID  model. 
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RAOIUS  taut  » 


Fig.  4.  Scalier  plot  or  temperature  versus  radius  Tor  the  same  prob- 
lem as  shown  in  Fig.  (2)  and  (3).  This  calculation  is  performed  with 
a grid  spacing,  in  both  r and  Z,  which  is  one  half  the  previous  value 
(twice  as  well  resolved). 
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LOG  (TIME) 


Ei*.  6 Central  temperature  as  a function  of  time.  The  initial  conditions  for  this  calcula- 
tion are  similar  to  those  shown  in  Fig.  (2)  except  A*,,  — 0.2  cm  and  R„  — 0.067  cm.  In 
the  previous  calculation.  Fig.  (S),  the  exterior  region  did  not  burn.  This  calculation, 
however,  includes  reactions,  conductivity  and  diffusion.  Thus  a flame  front  propagates 
from  the  interior  to  the  outside  wall.  This  front  is  coincident  with  the  location  of  the 
maximum  fluid  divergence. 
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Rig.  7.  Compulation  time  required  to  calculate  the 
diffusion  velocities  Tor  all  species  at  all  grid  points  (with 
the  number  of  grid  points  shown  along  each  curve).  As 
can  be  seen,  the  compulation  time  increases  more  slowly 
than  the  number  of  grid  points  as  vector  overhead  di- 
minishes in  importance. 


